package fem2;

import java.util.ArrayList;

import junit.framework.Assert;
import math2.delaunay.AbstractDelaunay;
import math2.delaunay.Delaunay2;

import org.junit.Test;

import fem2.enu.MeshPartType;

public class XfemIntegrationTest {

	@Test
	public void test1() {
		int p = 5;
		IntegrationManager im = IntegrationManager.getInstance();
		MeshPartType mpType = MeshPartType.QUAD4;
		double I = 0.0;
		int ng = im.getNumIntegrationPoints(mpType, p);
		for (int k = 0; k < ng; k++) {
			double[] xi = im.getIntegrationPoint(mpType, p, k);
			double w = im.getIntegrationWeight(mpType, p, k);
			I += w * Math.exp(xi[0] + xi[1]);
		}
		Debugger.watch(I);

		double exact = Math.pow(Math.exp(1) - Math.exp(-1), 2);
		Assert.assertEquals(exact, I, 1e-6);
	}

	@Test
	public void test2() {
		AbstractDelaunay d = new Delaunay2(-2, 2, -2, 2);

		d.addPoint(-1, -1);
		d.addPoint(1, -1);
		d.addPoint(1, 1);
		d.addPoint(-1, 1);
		d.addPoint(-1, -0.2);
		d.addPoint(1, 0.2);
		d.addPoint(0.0, 0.0);

		int p = 5;
		ArrayList<IntegrationPoint> ips = IntegrationManager.getInstance().createIntegrationPoints(
				p, d.getTriangles());

		double I = 0.0;
		for (IntegrationPoint ip : ips) {
			double[] xi = ip.getPoint();
			double w = ip.getWeight();
			I += w * Math.exp(xi[0] + xi[1]);
		}
		Debugger.watch(I);

		double exact = Math.pow(Math.exp(1) - Math.exp(-1), 2);
		Assert.assertEquals(exact, I, 1e-6);
	}

	@Test
	public void test3() {
		AbstractDelaunay d = new Delaunay2(-2, 2, -2, 2);

		d.addPoint(-1, -1);
		d.addPoint(1, -1);
		d.addPoint(1, 1);
		d.addPoint(-1, 1);
		d.addPoint(-1, -0.2);
		d.addPoint(1, 0.2);
		d.addPoint(0.0, 0.0);

		// double exact = 4.9999441603522285504368483088910579681396484375;

		Debugger.watch("using sub-triangles");
		for (int p = 1; p < 6; p++) {
			ArrayList<IntegrationPoint> ips = IntegrationManager.getInstance()
					.createIntegrationPoints(p, d.getTriangles());

			double I = 0.0;
			for (IntegrationPoint ip : ips) {
				double[] xi = ip.getPoint();
				double w = ip.getWeight();
				double r = Math.sqrt(Math.pow(xi[0], 2) + Math.pow(xi[1], 2));
				I += w * 1 / Math.sqrt(r);
			}
			Debugger.watch(I);
		}

		Debugger.watch("using normal quadrature");
		for (int p = 1; p < 9; p++) {
			MeshPartType type = MeshPartType.QUAD4;

			int ng = IntegrationManager.getInstance().getNumIntegrationPoints(type, p);
			double I = 0.0;
			for (int j = 0; j < ng; j++) {
				double[] xi = IntegrationManager.getInstance().getIntegrationPoint(type, p, j);
				double w = IntegrationManager.getInstance().getIntegrationWeight(type, p, j);
				double r = Math.sqrt(Math.pow(xi[0], 2) + Math.pow(xi[1], 2));
				I += w * 1 / Math.sqrt(r);
			}
			Debugger.watch(I);
		}
	}
}
